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Abstract: 

A minimal ofF-lattice model for a-helical proteins is presented. It is based on hydrophobicity forces 
and sequence independent local interactions. The latter are chosen so as to favor the formation of 
a-helical structure. They model chirality and a-helical hydrogen bonding. The global structures 
resulting from the competition between these forces are studied by means of an efficient Monte Carlo 
method. The model is tested on two sequences of length N = 21 and 33 which are intended to form 
2- and 3-helix bundles, respectively. The local structure of our model proteins is compared to that 
of real a-helical proteins, and is found to be very similar. The two sequences display the desired 
numbers of helices in the folded phase. Only a few different relative orientations of the helices are 
thermodynamically allowed. Our ability to investigate the thermodynamics relies heavily upon the 
efficiency of the used algorithm, simulated tempering; in this Monte Carlo approach, the tempera- 
ture becomes a fluctuating variable, enabling the crossing of free-energy barriers. 

Key words: protein folding, hydrophobicity, simulated tempering, global optimization, helical bun- 
dle, alpha-helix 
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1 Introduction 



Protein folding may be described on all levels of complexity, ranging from most simplistic concepts 
that do not even incorporate the concept of geometry, to all-atom representations including solvent. 
Evidently, different levels of descriptive complexity address different aspects of the protein folding 
problem. 

Addressing the thermodynamics of a given protein model is generally very difficult due to the free 
energy barriers present for compact chains Therefore, past thermodynamic studies of proteins 
have mainly been performed using lattice models. Lattice models lend themselves to thermodynamic 
calculations since the conformation space is discrete, enabling the enumeration of all conformations 
and therefore the calculation of exact results for short chains. For example, the minimal HP model 
of Lau and Dill |^ has been examined in quite some detail However, the approximations 

involved in lattice models are far from well understood |^ , and their geometry is clearly not that of 
real proteins. Furthermore, the energy barriers between the states may be poorly represented due 
to the discreteness of the conformation space. 

In this paper, we study the thermodynamic behavior of a simple off-lattice model for protein folding. 
Studies of similar, minimal, off-lattice models have been performed before, see e.g. 0, |[ ^, ^ 
The main difference between our model and those studied previously is that it is deliberately 
constructed so as to model the Ca backbone geometry of a-helical proteins. A somewhat similar 
approach was taken in [|l2|. However, these authors used a conventional Monte Carlo method, and 
reliable thermodynamic averages for global quantities were not obtained. 

Our starting point is the model recently suggested by Irback et al. [jll]. The major change is a 
modification of the local interactions which we here choose so as to model a-helical structure. We 
stick to the original concept of presenting each amino acid as a single site. Furthermore, only two 
types of residues are considered, hydrophobic and hydrophilic. The residues are linked by rigid 
bonds. The local interactions in our model are meant to model the chirality of the amino acids and 
the hydrogen bonds found in real a-helices. 

It has been shown by Kamtekar et al. that the binary sequence pattern of hydrophobic and 
hydrophilic amino acids is of central importance in the design of de novo proteins. The model pre- 
sented here has some potential to be of help when pursuing this type of protein design. In particular, 
it may give valuable hints on different thcrmodynamically stable structures of a given hydropho- 
bic/hydrophilic sequence; such information would certainly be useful in this design approach. 



2 Methods 



2.1 The Model 

Our starting point is the model proposed in to be referred to as the AB model, which is 
briefly described in the appendix. In the AB model, there are sequences which are thermodynam- 
ically stable at kinetically acceptable temperatures Also, there are strong regularities in the 
local structure of the chains, qualitatively similar to those for real proteins. However, at a more 
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quantitative level, the local structure observed in ||ll|] is very different from that of real proteins 
- partly because it is insensitive to space reflections (no chirality). The modifications of the AB 
model presented in this paper are meant to model a-helical local structure; in particular, chirality 
is introduced. Globally, the modifications lead to sizes similar to those of real proteins. 

Let us call the right-handed a-helical structure found in real proteins an ideal a-helix; this helix 
has 3.6 amino acids per turn and a translation of 5.4A per turn. An ideal a-helix is completely 
described by the geometry of three successive (virtual) Cq-Cq, bonds; that is by one torsional (a^) 
and two bend angles (Ti,ri+i). Equivalently, one may describe an ideal a-helix by the mutual 
distances of the four participating Ca atoms (plus the correct chirality). 

In our model, each residue is represented by a single site corresponding to the Ca position in the 
polypeptide backbone. These sites are linked by rigid bonds, bi, of length 3.8A. The shape and 
energy of an 7V-mer is specified by the TV — 1 bond vectors hi . Vij denotes the distance between 
residues i and j. (Ji, . . . ,aM is a binary string that specifies the primary sequence; we consider 
hydrophobic ((7 = 0) and hydrophilic ((7 = 1) monomers only. The energy function is defined to be 

The first term in this equation consists of Lennard-Jones (LJ) interactions. The LJ parameters Aij 
are chosen as A,j = 2-^/^r°-j for I J— 1 1< 4, where denotes the distance in an ideal a-helix. This 
choice of Aij gives the corresponding terms in eq. |l| a minimum at rfj. For all other (i, j), we set 

Aij 6.OIA, which is calculated from the average volume, 16lA^, of an amino acid Hence there 
are four different Aij values, three of which are determined by the geometry of an ideal a-helix, the 
fourth being defined by the average size of an amino acid. Finally, the depth of the LJ potential is 
chosen so as to favor the formation of a hydrophobic core; e(0, 0) = 1 for a hydrophobic-hydrophobic 
pair, e=l/2 for all other pairs. 

The sequence-independent local interaction Ei{i) depends on the torsional and bend angles of three 
successive bond vectors (a^, r^, r^+i), and is given by a negative gaussian with unit depth and width 

i^.(.)^^exp(- (-'--^)^ + (--^;f + (--+^'^)^ ) i.,(.)ehl,0] (2) 

where a and f are the torsional and bend angles in an ideal a-helix. Ei(i) is close to minus one if 
and only if the three participating bond vectors are close to a-helical structure, which is illustrated 
in Fig. |l]. The width is taken to be u; = 1/4 radw 14.32°; this choice is justified by the resulting angle 
distributions, which will be presented in section 3.1. The aim of the local interaction is twofold. 
First, it introduces a right-handed chirality into the model that is not present in the LJ term alone. 
Second, it is meant to model the hydrogen bonds of the a-helical conformation. Observe that the 
depth of the local interaction Eiii) equals the depth of a hydrophobic-hydrophobic contact. Exact 
values of the parameters and technicalities, as well as a comparison to the AB model, can be found 
in the appendix. 

In section 3.1, it will be shown that, not surprisingly, the model indeed gives rise to a-helical struc- 
ture. This is done by comparing the local structure of our model to that of real a-helical proteins. 
Quantitatively very similar results are found. Furthermore, we examine the overall topology of the 
simulated sequences. As will be shown in section 3.2, we find topologies corresponding to 2- and 
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Figure 1: Local energy Ei{i) (Eq. ^ displayed in one bond and one torsional angle. The second 
bend angle is set to f , the bend angle of an ideal a- helix. Ei(i) has a single minimum of depth one 
around the a-helical conformation and is close to zero elsewhere. Equipotential lines for Ei(i) e 
{—0.1, —0.2, —0.9} are displayed with circles. 

3-helix bundles. By looking at the Ca-Ca distance distribution functions, it can be seen that even 
the linear dimensions of our model proteins are similar to those of real ones. 



2.2 Sequences 

We have studied the behavior of two sequences in this model. These were deliberately chosen so as 
to be consistent with a-helical structure. Specifically, we start from a sequence segment of length 
14, which is known to occur in amphiphilic a- helices p5| . This segment has been used earlier for the 
design of (real) de novo a-helical proteins [|l3j. We take the first 9 positions of this A^=14 segment 
as our basic building block. Our two sequences are then obtained by taking two (three) copies of 
this building block, and connecting them by hydrophilic segments of length three. These sequences 
have length iV = 21 and 33 respectively, and are given in table 0, where the spaces are meant to 
clarify the construction of the sequences. It should be obvious that these are meant as templates 
for 2- and 3-helix bundles, respectively. 



Length N 


sequence ai, ...,aN 


21 


101100110 111 101100110 


33 


101100110 111 101100110 111 101100110 



Table 1: The two sequences studied. Basic building blocks of length 9 are connected by hydrophilic 
segments of length 3. 



2.3 Monte Carlo Methods: Simulated Tempering 

We calculate thermodynamic properties by using the method of simulated tempering [|l^, [l^, p^ . 
To simulate proteins is notoriously difficult, due to the presence of a rugged energy landscape. One 
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tries to overcome this problem in simulated tempering by treating the temperature as a fluctuating 
variable. This means that one simulates a joint distribution in conformation and temperature, which 
is taken to be 

Pib, k) cc exp(-5fc - E(b, a)/Tk) , (3) 

where Tk, k ~ 1, . . . ,K, are the allowed values of the temperature. The gk^s are free parameters 
that determine the probabilities P{Tk) of visiting the different temperatures. Although the method 
is free from systematic errors for any set of gk, it is crucial for the performance of the algorithm 
to make a careful choice of these parameters. In our simulations, the parameters gk were adjusted 
so as to have a roughly uniform distribution P{Tk). This was done by means of trial runs. Details 
may be found in [ p^ . 

The simulations reported in this paper were carried out using a set K — 20 temperatures, ranging 
from Ti = 0.15 to Tzo = 1.0 with l/Ti - l/T,+i constant. On a DEC ALPHA 200 (266 MHz), 
they took around 30 hours for = 21 and around 500 hours for A^ = 33. In both cases, around 
20 percent of the total computing time was spent on tuning of the weights gk ■ We also simulated 
a sequence of length A^ = 45, whose construction is analogous to that of the other two sequences. 
We spent 1000 CPU hours on this sequence, and believe that we obtained reliable estimates for its 
local properties. However, the results for global properties must be interpreted with care for this 
sequence. By contrast, we feel very confident that the simulations for A^ = 21 and A^ = 33 are under 
control. 



3 Results 



The overall thermodynamic behavior of our model turns out to closely resemble that of the AB 
model ||ll|. In particular, the chains compactify gradually with decreasing temperature, and the 
folding transition takes place in the compact phase. Furthermore, we do not observe bimodal 
distributions in energy, compactness or the local interactions Ei(i). All these observations are in 
agreement with earlier results for the AB model. By contrast, our model differs significantly from 
the AB model when it comes to structural properties. We shall therefore focus on these. Local 
structure as well as global topology and stability will be discussed in some detail. 



3.1 Local Structure: Comparison with Real Proteins 

In this section we examine the local structure of the chains in our model, and compare it to that 
of real a-helical proteins. The protein data were extracted using the Ca coordinates for 27 proteins 
containing mainly a-helices and almost no /3-sheets. All these 27 proteins are listed in the all- 
a-helical class of the SCOP database |lq|. The corresponding structures were taken from the 
Brookhaven Protein Data Bank (PDB) |l9| n. One of these proteins, the A^ = 54 DNA-binding 
protein lenh, which contains three a-helices [^0|, will be considered in more detail. 

The local structure will be probed in two ways, by bend and torsional angle distributions and by 
bond-bond correlations. 

^The PDB accession codes are: 3sdh, Ictj, lenh, 2erl, 2end, His, Ihme, Infn, Ibcf, Irhg, lacp, Irop, Icoo, 4icb, 
lutg, Ifia, 2wrp, 2tct, Ifps, lecm, laep, laxn, Ihiw, 2abk, leci, 2abd, lc5a. 
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Figure 2: Bond (r^) and torsional {ai) angle distributions. Left: For the 27 a-helical proteins taken 
from the PDB database. Right: Model results for the iV = 21 sequence at r = 0.15. The other two 
sequences simulated, with iV = 33 and N = 45, give very similar distributions. Both scatter plots 
contain 1000 data points. 



In Fig. |2[ left plot, bend and torsional angle distributions are displayed for the set of 27 a-helical 
proteins. The dominance of a-helical structure manifests itself in the pronounced peak centered at 
a « —129° and r « 90°. The same distributions, for our iV = 21 sequence, are plotted on the right 
side of Fig. The results for the other two simulated sequences are very similar. This is in line 
with the behavior of the AB model, where the angle distributions are fairly sequence-independent 



as well |11 . Not surprisingly, we find that there is indeed a peak in the distribution centered at 
the position of the ideal a-helix. This peak contains roughly the same amount of probability as for 
real a-helical proteins. Moreover, the width of this peak is very similar to that of the experimental 
structure, which justifies our ad hoc value for the width w of the exponential in Eq. |[ Another and 
much weaker peak in the angle distribution for our model is found at (a, t)=(0°, 90°), which is an 
artifact of our model not present in real proteins. 

Our second way to monitor local structure arc the bond-bond correlations, which are studied using 
the function 

^"(^^^^ nr^'i ) (4) 
\\h\\h+d\/ 

where (•) denotes the average over all positions i and observed structures (with d fixed). Cb{d) is a 
measure of the average alignment of bond vectors at a given topological distance d along the chain. 
It is normalized so that Cb(0) = 1. In Fig. |3| (dotted line, left plot), we show the correlation function 
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Figure 3: Left: Correlation functions Cb{d) for our model at T = 0.15 for N — 21 (solid line) and 
A^ = 33 (dashed line). Cb{d) for the 27 a- helical proteins listed in the text is given by the dotted 
line. Right: Cq-Cq, distance distribution P{r) for the (molecular dynamics refined) experimental 
structure of the DNA binding protein lenh (TV = 54, solid line), and our model (TV = 45, dashed line) 
at T = 0.15. Nearest-neighbor pairs are not included. 

Cb{d) for the 27 a- helical proteins. As can be seen from this figure, there are significant correlations 
at least out to separations of about ten residues. The oscillations can be related to the presence of 
a-helical structure, which has a period of 3.6. In Jill it was shown that the AB model gives roughly 
the right correlation length, but fails to reproduce the periodicity of 3.6. Cb{d) for our model is 
shown in Fig. ^, using A^ = 21 (solid line) and N = 33 (dashed line). Although the correlations decay 
somewhat faster in the model, the results are very similar to those for the real proteins. 



3.2 Global Structure and Stability 

It is an essential feature of any plausible model for protein folding that it allows for structural 
stability. In the following section, we examine our model with respect to this important aspect 
using two methods. First, we study structural stability by measuring mean square fluctuations, 
(5^). Second, the structure of the folded conformations will be characterized by using a measure for 
a-helix formation. 

A commmonly used measure of the similarity between two conformations a and b is the mean square 
distance 5^^, which is defined as 



where \xl ~ x\ \ denotes the distance between the sites and xl , and where the minimum 
is taken over translations and rotations. The probability distribution P{5^), measured on a ther- 
modynamic ensemble with fixed sequence and fixed temperature, gives valuable information about 




(5) 
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Figure 4: The function Ei (i) of Eq. |^ for the experimental structure of DNA binding protein, 
PDB accession code lenh (A^ = 54). One easily identifies three a-helices from this plot; for a-helical 
structure Ei{i) is close to minus one. The resulting a-helix assignments are in agreement with those 
in the PDB entry. 



structural stability . In particular, if there is a number of relatively well-defined structures present, 
the distribution will have a narrow peak close to w 0; in this case the ensemble contains many 
pairs of similar conformations. 

Our second method for characterizing structure and stability is to monitor the formation of a-helical 
structure. A convenient way of doing this is to utilize the function Ei{i) (Eq. ||); besides using it for 
modelling hydrogen-bonding in the a-helical structure, Ei(i) can also be used to detect a-helices 
in a given structure. To illustrate this, we plot Ei{i) for the experimental structure of the DNA 
binding protein lenh in Fig. ^ From this plot, it is immediately clear that the structure contains 
three a-helices. The first a-helix is identified by the drop of Ei{i) from zero to slightly above minus 
one in the region from position i = 7 to i = 19. It should be observed that 16 amino acids participate 
in this helix; the positions of four successive C^-atoms determine one value of Eiii). Below, we will 
perform a similar analysis for our model proteins. 

For the A^=:21 sequence at T=:0.15, the distribution P{6^) has a mean value of ((52) = 3.6±0.5A . 
This result may be compared to those of ||l^ for the AB model. There, six different A^= 20 sequences 
were studied, and (5^) varied between lA and lOA . It should be pointed out, however, that the 
size of the monomers, as measured by the dimensionless ratio Aij/r^^i+i, is slightly larger in the 
present model (it is equal to one in the AB model). 

One could be tempted to conclude from the P{5'^) distribution, which is shown in Fig. || (right, 
dotted line), that the A^ = 21 system has two similar but distinct states. A closer analysis reveals 
that there are actually four distinct states; the system spends almost all time, 95%, in the near 
vincinities of these. This can be shown by measuring the simultaneous 6^ distances to all these four 
states. The mutual distances between these states range from l.sA to 3.7A . When it comes to 
helix content, these four states are very similar. The function Ei(i), as displayed in Fig. ^, makes it 
evident that the overall structure is that of a two helix bundle. We also looked at the expectation 
values of the bend and torsional angles and their variances, which showed that the helices are indeed 
very stable and common to the four different states. The structure of one of the four states is shown 
in Fig. H (the other three are very similar). From Fig. ^ it can also be seen that (1) we succeed to 
obtain a hydrophobic core and (2) the a-helices show the patterning intended in the construction 
of the sequences. 
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Figure 5: Representative structures for the most probable states for iV = 21 (left) and = 33 
(middle). Hydrophobic and hydrophilic amino acids are represented by black and white spheres, 
respectively On the right side, P{S^) is plotted for (1) N = 21, r = 0.15 (dotted line), (2) iV = 33, 
r = 0.15 (solid line), and (3) iV = 33, r«0.40 (dashed line). 




Figure 6: Thermodynamic averages of the local interactions Ei{i) a,t T — 0.15 for A^ = 21 (left) 
and A^ = 33 (right). These curves clearly show that the dominating states for iV = 21 and A^ = 33 
contain 2 and 3 helices, respectively. This is also supported by the fact that fluctuations of bend 
and torsional angles are small in the helical regions. 

The P{S'^) distribution for = 33, see Fig. ^ (right plot), also has a pronounced peak at w 0. 
However, a significant amount of probability is also found around « 15A . The existence of this 
outlier in P((5^) shows that the overall geometry must be different for the different states present. A 
detailed analysis, similar to that for A^=:21, shows that P{S^) is dominated by three states having 
mutual distances of around 13.8, 15.0 and 16. 3A . Taken together, these three states contain 97% 
probability at T = 0.15, distributed as 68%, 17% and 12%. These three states share three common 
a-helices; the positions of the helices along the chain can be seen from Fig. ^. How do the three 
states differ? For the most probable state (68%), the three helices are coordinated in such a way 
that the overall topology has right-handed chirality , as can be seen from Fig. ||, middle plot. For 
this state, 26 monomers participate in the three helices. The other two states are described as 
follows: the first two helices and the connecting loop form a U-shaped entity (the helices are not 
really parallel). The third helix points back into the U, for one state on top and for the other below 
the U-plane. The right-handed state is shghtly less probable (12%) than the left-handed one (17%). 
The helices of the most probable state are present in these two states as well; however, the helices 
of these two states are longer with 31 and 33 monomers participating in the sense that 22 and 24 
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values of Ei{i) are close to minus one. 

Knowing about this three-fold structural degeneracy could be of interest forprotein design. If the 
design is done in terms of the hydrophobicity pattern, as for example in ||l3|, the designer could 
easily knock out two unwanted states by introducing unfavorable residues. 

As stated in section 2.3, we also tried to simulate a sequence with N = 45, intended to design 
4-helix topology. 1000 CPU hours were insufficient to obtain fully reliable results. However, we 
feel confident that we managed to adjust the weights gk pretty well and obtained a selection of 
low-energy structures that we rate as representative when it comes to local structure and overall 
size of the chain. To give an idea about the size of our model chains, we display the C'a-Ca distance 
distribution function in Fig. |^ (right figure, dashed line). In this graph, we compare it with the 
same distribution for the DNA binding protein lenh, which contains 54 amino acids. The general 
shape of the two curves is similar, taking into account the shorter length of the model protein. This 
shows that even the size of our model proteins resembles that of real ones. 



4 Conclusions 



A minimal off-lattice model for a-helical proteins has been presented, a-helical structure is obtained 
by setting local LJ parameters to values consistent with a-helical geometry and introducing local 
interactions modelling chirality and a-helical hydrogen bonding. The global LJ parameters are 
determined by the size of an average amino acid. 

Thermodynamic results for two sequences of length 21 and 33 have been presented. The construction 
of these sequences was inspired by the method of Kamtekar [|l3j to design real de novo proteins. This 
turns out to be a successful strategy here as well, since we observe a-helices at the intended positions 
along the chains. Furthermore, the helices observed in our model show the intended hydrophobicity 
pattern with one side being hydrophobic and the other hydrophilic; the hydrophobic sides of the 
helices point inwards. 

The local structure has been explored using angle distributions and bond-bond correlations. When 
we compare results from our model with those for real a-helical proteins, these observables are very 
similar. 

By inspecting the mean square distribution, P{5'^), it has convincingly be shown that the A'' = 21 
sequence is structurally stable. The overall topology was found to be that of a 2-helix bundle. For 
A^ = 33, we found three globally different states which dominate the low temperature phase. These 
three states have three a-helices in common and differ from each other by the orientation of these 
hehces. 

Simulated tempering is the key to our ability investigate the low-temperature thermodynamics. The 
simulations took 30 CPU hours for A^ = 21 and 500 CPU hours for A^ = 33. 1000 CPU hours were 
not sufficient to obtain reliable results for A = 45, but we feel confident that we managed to properly 
adjust the weights gk- Therefore, we believe that simulating the A^ = 45 sequence in a reliable way 
will be possible in the near future. Also, we feel confident that the simulations can be speeded up 
by further algorithmic improvements, such as a better way of distributing the temperatures Tk (see 
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Figure 7: Illustration of an ideal a-helix: The polypeptide backbone has 3.6 amino acids per turn 
with a pitch of 5.4A, as indicated by the height of the cylinder. The cylinder has a diameter of 
4.5. .A. 



e.g. m 



5 Appendix 

5.1 The Ideal a-Helix: Model Parameters 

Five parameters of our model are determined by the geometry of an ideal a-helix. A sixth parameter 
is determined by the average size of an amino acid. In this section, the calculation of these parameters 
is presented. 

A discrete helix is uniquely described by its bend and torsional angles, or, equivalently, by its 
pitch, the number of entities per turn and its chirality. An ideal a-helix has a translation of 5.4A 
and 3.6 amino acids per turn with right-handed chirality, as is illustrated in Fig. ^ This is easily 
reformulated in terms of bend and torsional angles which we denote by f and a, respectively. 

In order to make everything well-defined, we have to give an exact definition of bend and torsional 
angles. Three succesive bond vectors hi, bi+2 define the two bend angles Ti^n+i and the torsional 
angle a^. The bend angle r is defined to be 180° if the two involved bond vectors are parallel. The 
torsional angle a is taken to be the angle between the two vectors bi x bi^i and x 6^+2. Therefore, 
Q! = if and only if the three involved bond vectors are situated in a plane, a ranges from —180° to 
180°, and is chosen to be negative for a right-handed and positive for a left-handed system. 



A single LJ term in Eq. |l|. 




(6) 



has its minimum of depth e at = Aij ■ 2^/^; the minimum can be moved by changing Aij. We 
adjust Aij for | i — j |< 4 such that the minimum is situated at the distance observed in an ideal 
a-helix, r°j. For all other pairs, — i |> 4, we choose the LJ parameter Aij by tuning the volume 
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of the LJ sphere to that of an average amino acid, I6IA : 

|z-j|> 4 ^ Ay = 2-1/^ X {j^l&\k^ « 6.01A (7) 

To summarize, one obtains the following values for f, d, r^- and Ay-: 

f«90.53° a«-129.61° 

|j-zh2 |j-»h4 |j-z|>4 
[A] 3^8 5i)4~ 097 

Ay [A] - 4.80.. 4.49.. 5.52.. 6.01.. 



5.2 The AB model 

The a-helical model presented in this paper is derived from the 3D AB model suggested recently 
by Irback et al. jll), which in turn is closely related to the 2D model proposed by Stillinger et. 
al H, |l^. To make comparisons easier, we here give the energy functions for both the AB model 
and the present model; 



N / , . -.12 /4\6\ N-2 r r N-3 , , 



t+2\ 



2^2 



i=l j=i+2 



They differ only in the choice of Aij and in the local interactions; for the AB model, Aij = 3.8A 
independent of i and j. The dependence of the AB model on the strengths of the local interactions 
(ki,K2) was discussed to some extent in [|ll|, where it was found that the local interactions are 
necessary in order to obtain regularities in the local structure and thermodynamic stable sequences. 
The final choice of (ki,K2) in |ll[ was {ki,K2) = (—1,0.5). 
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